Electronic topological transition in sliding bilayer graphene 
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We demonstrate theoretically that the topology of energy bands and Fermi surface in bilayer 
graphene undergoes a very sensitive transition when an extremely tiny lateral interlayer shift occurs 
in arbitrary directions. The phenomenon originates from a generation of an effective non-Abelian 
vector potential in the Dirac Hamiltonian by the sliding motions. The characteristics of the transition 
such as pair annihilations of massless Dirac fermions are dictated by the sliding direction owing to 
a unique interplay between the effective non-Abelian gauge fields and Berry's phases associated 
with massless electrons. The transition manifests itself in various measurable quantities such as 
anomalous density of states, minimal conductivity, and distinct Landau level spectrum. 
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I. INTRODUCTION 

Changes in the topology of Fermi surfaces known as 
Lifshitz transition^ alter physical properties of metals 
signifcantly^"'^ . Though such an electronic topological 
transition (ETT) has been pursued in various materials, 
its realization requires large external perturbations such 
as alloying or applying high pressure that hinder clear de- 
tections of the transition Recent progress in measur- 
ing low energy electronic structures of bilayer graphene 
(BLG)^"^ provides a new opportunity to explore the ETT 
because of its unique electronic structures and because of 
the possible noninvasive control of chemical potential of 
the system^. 

In BLG, two coupled hexagonal lattices of carbon 
atoms are arranged according to Bernal stackingiSrJ^. 
Because BLG has a large degeneracy at the charge neu- 
tral poiniiiSr— , there have been intense discussions on 
possible many-body effects in the system^i^iii^"— More- 
over, since electrons in a single layer graphene (SLG) 
behave as relativistic massless fermions^^, BLG provides 
a unique playground to control interactions between rel- 
ativistic particles coupled with the relative mechanical 
motions of two layers. Hence, the effects of rotational 
stacking fault on physical properties of BLG have been 
studied extensively^i^"— : however, the effect of sliding 
one layer with respect to the other has not. This mechan- 
ical motion is important because the interactions between 
the two layers are sensitive to the deviation from Bernal 
stacking22r— and extremely small sliding will change its 
low energy electronic structures significantly. 

In this paper, we predict a very sensitive topological 
change in the energy bands and Fermi surfaces of BLG 
when sliding motion or interlayer shear occurs. It is 
demonstrated that a peculiar coupling between the ef- 
fective gauge potential with SU(4) symmetry generated 
by sliding motions and Berry's phase of massless Dirac 
fermions play a crucial role to change the topology of 
low energy bands of BLG. It gives rise to either pair an- 
nihilations of massless Dirac fermions or generations of 



fermions by absorbing fermions with topological charges, 
depending on sliding directions. This will offer new op- 
portunities to realize the ETT driven by non-Abelian 
gauge fields with gentle maneuverable mechanical oper- 
ations. 

We start with a detailed description of our first- 
principles calculation methods including a correction for 
interlayer dispersive forces. The energetics and changes 
of the interlayer distance for sliding BLG are presented 
also. Then, the low energy electronic structures obtained 
by the calculation are discussed when BLG experiences 
very small sliding between the two layers. The next three 
sections introduce a model Hamiltonian for the system 
and discuss the role of gauge potential and Berry's phase 
for the changes in low energy electronic structures. Sev- 
eral spectroscopic consequences will be discussed in the 
final section. 



II. FIRST-PRINCIPLES CALCULATION 
METHODS AND ATOMIC STRUCTURES 

Our electronic structure calculation employs the first- 
principles self-consistent pseudopotential method^'^ us- 
ing the generalized gradient approximation (GGA) for 
exchange-correlation functional emplemented by Perdew, 
Burke, and Ernzerhof (PBE)^^. A kinetic energy cutoff 
for wavefunction of 65 Rydberg is employed and a plane- 
wave basis set is used. The ion core of carbon atom 
is described by ultrasoft pseudopotential^^. A fc-point 
sampling of 60 x 60 x 1 A: points uniformly distributed in 
the two-dimensional Brilouin zone (BZ) is used in self- 
consistent calculations and 150 x 150 x 1 k points is sam- 
pled to obtain electronic energy bands on a rectangular 
grid of 1.2% of the first BZ size around the if-point. 
Since we have dealt with sub- Angstrom displacements of 
atoms, we have tested our calculation as increasing the 
cutoff to 120 Rydberg (corresponding to a kinetic energy 
cutoff of 480 Rydberg for charge density and potential) 
finding no differences in results. The total charge was 
calculated by using the Marzari-Vanderbilt cold smear- 
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FIG. 1: (color online) (a) In the top (black) and bottom layer 
(green), two sublattices are denoted by empty (^i) and filled 
(Bi) circles and by empty (B-i) and filled (^2) rectangles, 
respectively. The bottom layer (layer 2) slides with repect 
to the top (layer 1) by ds- In the top, 5i = (0, — l)ac, &2 — 
(-^,i)ac and ^3 



2-, jy'^c ".iiu — (-^, |)ac. (b) Schematic diagrams for 
stacking under various sliding directions, (c) The nn intra- 
and inter-layer hopping parameter, —70 and —71. (d) The 
nnn inter-layer hopping (between Bx-A'i), —73 with sliding, 
(e) The smallest inter-layer hopping (between A1-A2 or Bi- 
B2), —74 with sliding. 



ing schemers. All atomic coordinates are relaxed and the 
nearest-neighbor carbon-carbon distance (oc) in a single 
layer of graphene (SLG) is found to be 1.425 A. We also 
perform first-principles calculations again for all differ- 
ent stacking geometries by using another computational 
package^l finding no difference. 

Since the GGA cannot describe the interlayer inter- 
action between graphene properly, we have employed a 
seimemprical addition of van der Waals (vdW) forces to 
our calculations following Grimme's proposal^. Within 
the semiemprical method, the total energy of the system 
is i?pBE-i-vdw — EpBE + Evdw where -Epbe is the total 
energy from GGA functional of PBE and the total energy 
given by dispersive forces is Ey^yf which can be written 

as -BvdW = -5 Csjj J2r /dampd^y + R\)\nj + R\^^ , 

where /damp = se • (H-e"''-(l-^-+«l/'-«-i))^i, f„- = f,-r, 
is a vector for carbon-carbon distance, R a lattice vec- 
tor, sg a scaling parameter, and d a damping paramter, 
respectively^^i^. The coefficient of dispersive forces of 
Ceij and a sum of vdW radii tq are computed for each 
pair of atoms i and j such that Cgij — \J C^i Cqj and 
''o = roi+roj. We have used Cq = 1.75 J nm^ mol~^ and 
ro = 1.451 A for carbon atom suggested by Grimme^ 
and d = 20, sg = 0.65 for our GGA calculations. Here, 
i and j for -Evdw run through all atoms in the unit cell 
and R satisfies a criteria of \ fij + R\ < lOOA. 

Atomic structures of sliding BLG are shown in Figs. 
1(a) and 1(b). Without sliding, carbon atoms in BLG 
are arranged according to Bernal stacking - carbon atoms 
in one sublattice of the uppler layer are right on top of 
ones in the other sublattice of the lower layer. When 
the bottom layer slides with respect to the top layer, the 
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FIG. 2: (color online) Interaction energy per atom (Ei,) be- 
tween two graphene layers (a) with Bernal (AB) stacking and 
(b) with AA stacking as a function of interlayer distance 
(Ant). -EpBE, Evdw and -BpBE+vdw denote the interaction 
energies calculated by PBE functional only, vdW correction 
only and including both. Calculations results of i?pBE and 
i5pBE+vdw (retangles) are fitted to a second-order Birch equa- 
tion (solid lines). Eb = Etot — 2etot where Etot is the total 
energy of bilayer graphene per atom and etot that of single 
layer graphene per atom. In (c) and (d), the variations of 
binding energy (AE) and interlayer distance (ADint) with 
respect to those of AB-stacked bilayer are shown as functions 
of sliding distance in units of Oc along Si (red) and 5^ (blue). 
The negative AE denotes the decrease in binding energy. We 
set the binding energy and interlayer distance of AB-stacked 
bilayer to be zero in (a) and (b). 



sliding vector dg can be written as a linear combination 
of Si [i = 1, 2, 3) which connect the nearest neighbor (nn) 
carbon atoms in the top layer (For definitions of 5i {i — 
1, 2, 3), see Fig. 1). Considering the lattice structures of 
sliding bilayer, it is easy to check that the sliding along 
±(Si is equivalent to one along ±62(3)- We also note that 
a sliding along 83 is equivalent to one along —Si. Full 
sliding along Si brings an AB-stacking bilayer graphene 
to a BA-stacking bilayer graphene while one along 63 
shifts an AB-stacking bilayer to a AA-stacking bilayer. 

We find that the equilibrium interlayer distance of bi- 
layer graphene (BLG) in Bernal stacking is 3.348A and 
its binding energy is 38.76 meV (Fig. 2 (a)), which are 
in good agreement with other calculations and available 
experimental data^^Ti^. When bilayer graphene has AA- 
stacking (one layer is right on top of the other layer), 
the interlayer distance increases to 3.635A and bind- 
ing energy decreases to 30.77 meV (Fig. 2 (b)). When 
one graphene layer slides with respect to the other along 
either Si or S3 direction, the binding energy starts de- 
creasing and interlayer distance increases agreeing with 
a previous study^o (Fig. 2 (c) and (d)). We note that 
the change in the binding energy (< 0.19 meV) and in- 
terlayer distance (<0.01 A) is quite negligible when the 
sliding distance is less than 0.14 A (about 10 % of the 
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FIG. 3: (a) The small rectangle near around K-point in the 
first Brillouin zone is enlarged to show the three-fold sym- 
metric vectors koi = 1,0), kD2 = Pd(1/2, \/3/2), and 
kns = Pd(1/2, —y/3/2) with respect to koo = 0, where pn = 
(2/3ac)(7i 73/70)- (b) Energy contour for the valence bands of 
BLG without sliding. For all panels, = lOOx (fcx-FiT) /TK 
and Ky = 100 x ky/VK. The thick contour is an iso-energy 
{E — —1.1 meV) curve crossing the three saddle points. The 
blue (red) contours denote the hole (electron) pockets. Energy 
contours for the valence band of BLG with sliding ds — (c) 
0.02<5i, (d) 0T5i,je) -0.025i, (f) -O.Wi, (g) 0.012(^3 - (5i), 
and (h) 0.1(53 — Si). The contour interval is set by 0.5 meV 
for (a), by 1.0 meV for (c), (e) and (g) and by 10 meV for (d), 
(f) and (h). The thick contours that cross the saddle points 
are at E ^ (c) -7.2 meV, (d) -32.7 meV, (e) -10.8 meV, (f) 
-48.5 meV, (g) -9.7 meV and (h) -74.0 meV. 



bond length, Oc). 



III. LOW ENERGY BAND STRUCTURES 
FROM GGA CALCULATIONS 

Considering the low energy electronic structure of BLG 
in Bernal stacking, our calculations show, as in a previous 
study, four Dirac cones formed around the Fermi energy 
{Ep) at four Dirac points {koi , * = 0, 1, 2, 3) near around 



JC-point [Fig. 3(a) and 3(b)}iaiii. The magnitude of 
koi {i = 1,2,3) is about 0.4% of the distance between 
F- to A'-point (TK) [Fig. 3 (b)]. As the energy moves 
away from the Ep, the four Dirac cones merge to form 
three saddle points between the cones [Fig. 3(b)]. The 
calculation result can also be described by an effective 
HamiltonianiS, 

Hcff(fc) = hV3f-k+^{f* ■ k)T^{T* ■ k), (1) 
71 

where Va = (3/2)7Q,ac/fi. {a = 0,3), hk = h(kx,ky) is the 
crystal momentum from X-point, and r = (T^^Ty) and 
T* are Pauli spin matrices and their complex conjugates 
(See Fig. 1 for definition of 7^). Typical estimates are 
70-3 eV, 71 - 73 - 7o/10 and 74 ~ 70/2O (Ref^). 
We will neglect 74 and discuss its role in Sec. VL If we 
expand the effective Hamiltonian around each k^i, we 
have one isotropic Dirac Hamiltonian, JIdo — hv^r ■ Sko 
at kuo, an anisotropic Hdi = —fiV3{T*Skxi + 3T*Skyi) at 

A: £11, and two anisotropic others at k]j2 and kjj^ which 
can be obtained by rotating Hm by ±27r/3 respectively. 
Here, Ski = {5kxi,5kyi) = k - koi, [i = 0, 1, 2, 3). 

We find that the low energy bands of BLG change 
dramatically when one layer slides with respect to the 
other in an extremely small amount and in arbitrary di- 
rections. Let the bottom layer slide with respect to the 
top by ds{= {dx,dy)) [Fig. 1(a) and 1(b)]. First, when 
the bottom layer moves along Si direction by 0.028 A, i.e., 
ds = 0.02(^1, only two Dirac cones at k]j2 and k^^ remain, 
instead of four cones for BLG without sliding as shown 
in Fig. 3(c) The energetic position of the saddle points 
of the valence band decreases from —1.1 meV to —7.2 
mcV. When the bottom layer slides further {ds = O.lSi), 
the topology changes again and the saddle point energy 
significantly decreases to —32.7 meV [Fig. 3(d)]. Sec- 
ond, when the bottom layer moves along the direction 
opposite to the previous one {ds = — 0.02(5i), the low en- 
ergy bands changes again completely. In this case, the 
Dirac cone at fcui moves along the —k^ direction and an 
anomalous Dirac cone with sickle-shaped energy contours 
appears at a new Dirac point instead of three cones [Fig. 
3(e)]. If the bottom layer slides further by 0.14 A along 

— ^1 {ds = — 0.1(5i), the topology remains the same and 
the anomalous cone comes to have an anisotropic shape 
[Fig. 3(f)]. In this case, the saddle point energy decreases 
to -48.5 meV (almost 500% of pristine one) [Fig. 3(f)]. 
Finally, when we slide the bottom layer by 0.012((53 — Si) 
[Fig. 3(g)] and 0.1{S3 - Si) [Fig. 3(h)], the topological 
changes are similar to the case for the sliding along —Si. 
The saddle point energy for the sliding by 0.1(53 — Si) 
decreases dramatically down to —74.0 meV as shown in 
Fig. 3(h). After a comprehensive search for the topologi- 
cal changes by the sliding along arbitrary directions (not 
shown here), we find that the topology of saddle point en- 
ergy contours is all similar to those shown in Fig. 3(e)-(h) 
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(one crossing point) except the topologically distinctive 
phase in Fig. 3(c) (two crossing points). 



IV. EFFECTIVE MODEL HAMILTONIANS 
AND NON-ABELIAN GAUGE POTENTIAL 

The hypersensitive topological changes found by first- 
principles calculations in the previous section demon- 
strate that the sliding motion creates interactions be- 
tween the effective non-Abelian SU(4) gauge field back- 
ground and massless fermions. In the presence of a very 
tiny sliding, the nnn interlayer interaction (73) is not 
constant any more but depends exponentially on the dif- 
ferent pair distances between carbon atoms in top and 
bottom layers as shown in Fig. 1(d). The asymmet- 
ric inter-layer hopping interaction produces a constant 
pseudo-gauge potential in the terms of Hamiltonian con- 
taining 73 only, being similar with an effective Hamilto- 
nian for strained SLG^"— . Hence, the effective Hamil- 
tonian for sliding BLG can be written as, 

H',s{k) = hvsf- (fc - A) + ^(f* • fc)T,(T* • k) (2) 

Pd 

where = ^[v^/ {Hvq) and j[ is the reduced nn inter- 
layer interaction. (See Appendix for derivation of Eq. 
(2)). We find that the constant vector potential (A) ex- 
plicitly depends on the sliding vector, A = {Xx,Xy) = 
l3{dy, —dx)- By fitting the energy bands obtained by the 
model to ones by the first-principles calculation, we can 
estimate that /3 is about The gauge symmetry in 

Eq. (2) apparently seems to be broken since A is absent 
in the quadratic term of k. However, when we expand 
Eq. (1) at each Dirac point [koi], the four Dirac cones 
shift depending on both A and their own positions as 
km km + A,(\) {i = 0,1,2,3). If |A| <C pd, Ao = 



(Ax,Aj,), Ai — (— Aj,,-j-), A2 — {—^^^ 

and A3 = ( , ^•^^ — ^^). So, each Dirac cone moves 
along a different direction depending on its position. 

The essential feature of the system that enables the 
electronic topological transition in BLG is the shift of 
four Dirac points along different directions under the slid- 
ing of one layer against the other. As we have shown, the 
effect of sliding is to replace Ski — k—kjji with 6ki — Ai in 
the effective Hamiltonian for the low energy modes near 
the Dirac points and this feature suggests that the slid- 
ing in BLG induces a non-Abelian background gauge field 
associted with these modes. In order to show this more 
explicitly, we assume that before sliding, there are four 
massless Dirac fermions, labeled by i = 0, 1,2,3, whose 
Hamiltonians are all isotropic and of same chirality: 



5k., , 



0,...,3. 



(3) 



We can then form a quadruplet out of four fermions, 
so that the index i is now viewed as the "color" index 



of SU(4) symmetry, and combine the four Hamiltonians 
■Hui compactly asHo = hvz{f* ■ 5k) ® I, where I is the 
4x4 identity matrix. If we now introduce the background 
gauge field A — {Ax, Ay) for the SU(4) symmetry with 



/I 



A;^; = Xx 



^0 



-1 



Xy^ 

V3 



and 



A - — 



+ Xy 



1/ 



/I 



V 



so that 5kl — > Ski — A in the Hamiltonian, then Eq. ([3]) 
now becomes 



hvsT* ■ {5k, - Ai) , i = 0, . . . , 3, 



(4) 



where Ai are precisely the shifts of the Dirac points 
shown at the beginning of this section. Since both Aa, 
and Ay are linear combinations of the generators of 
SU(4), we now see the effect of sliding as if introducing 
a non-Abelian background gauge field. 

The low energy Hamiltonian of bilayer graphene is 
different from what we have just shown above in that 
Hdo has opposite chirality from TiDi for i — 1, 2, 3, and 
that the latter three Hamiltonians are anisotropic. A 
parity inversion Sky — > —Sky for "Hdo f^nd anisotropic 
rescalings of Ski for Hoi {i = li2,3) transform SU(4) 
isotropic quadruplet to the low energy Hamiltonian of 
BLG. This transformation is equivalent to multiplying 
the Pauli matrices Tx and Ty with (different) constants 
in Eq. @ — for example, {tx, Ty) — {tx, —Ty) for i = 0, 



and {Tx,Ty) 



-{Tx,STy) for i 



1 



thereby break- 



ing the SU(4) symmetry of the previous paragraph. But 
the structure of the Hamiltonians otherwise remains the 
same, and in particular, the shifts of the Dirac points are 
still given by A^ of Eq. dJ]). 



V. ROLES OF BERRY'S PHASE IN 
ELECTRONIC TOPOLOGICAL TRANSITION 

The characteristics of the ETT are ruled by unique in- 
terplay between the effective non-Abelian vector poten- 
tial and conservation of Berry's phase. The Berry's phase 
{(f>B) for each Dirac cone at km can be calculated by us- 
ing = §^dk-A{k) where A{k) = i{uDi{k)\dji\uDi{k)), 
r is a path enclosing each Dirac point, and U]ji{k) is 
a single-valued spinor-like eigenfunction of each Dirac 
Hamiltonian at kjji. Without sliding, we have (f>B = 
for the massless Dirac fermions around k^i {i = 1,2,3) 

and (pB = — TT for the ones at the center fc^o [Fig- 4(a)]. 
For a path enclosing all Dirac points at higher energy, we 
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FIG. 4: Three dimensional plots of low energy bands from the effective Hamiltonian of Eq. (1) with A = (a) Without 

sliding (A = 0), three anisotropic Dirac cones with Berry phase tt at koi (i = 1, 2, 3) denoted by dotted red arrows in all panels, 
and isotropic one with — tt at koo- Energy bands for sliding BLG with (b) A^; = — pd/5 (c) — pd/4, (d) —po/S, (e) Pd/5, (f) 
3pd/4 and (g) l.Spo. In term of sliding distance dy, (b) dy ~ -0.003 , (c) -0.004, (d) -0.005, (e) +0.003, (f) +0.012 and (g) 
+0.029 A. The solid red arrows in (b) and (e) indicate the movements of Dirac cones at sliding. 



can use the quadratic Hamiltonian (1) and 4>b is +27r— i^. 
So, the total </)b is always conserved to be 27r [Fig. 4(a))^. 

Now, to reveal the role of Berry's phase explicitly, let 
us consider exactly solvable cases without the assumption 
of |A| ^ pd- For sliding along ±6i direction, i.e., dg = 
(0,+dy) {dy > 0), the effective vector potential is given 
by A = (A„0) = iTpdy,0). Here, A, < (A, > 0) 
when sliding along +Si{—Si) direction. Then we have 
four different local vector potentials Ai for each kni such 
that 



An 




f (+l-,/l + l^,0|. 



A, 



We note that Aq + Ai ^ and A2 + A3 ^ 0. Thus, the 
Dirac cone at k^o and the one at fc^i move in opposite 
direction when bilayer graphene slides along ±Si and so 
do those at k£i2 and kjj^. It is also noticeable that koo + 
Aq — kui + Ai when A^; — — and koi + = 
ko2 + A2 — ko3 + A3 when A^; = 3pd/4. Hence, when 



the bottom layer slides along —y direction (i5i direction) 
by dy = pd/{4^I3), the two Dirac cones at kno and koi 
meet at (— p£)/2, 0) while three cones at k^i, ko2 and 
kD3 meet together at (pd/2,0) when sliding along +y 
direction {—61 direction) by dy = 3p£)/(4/3). 

For the sliding along +i5i direction, the effective Hamil- 
tonian shown in Eq. (2) can be expanded around fc+ = 
{~Pd /2, 0) so that the resulting Hamiltonian can be writ- 
ten as 



flV3Tx 



PD 



+ 2hv3TySky (5) 



where Sk ~ {Skx,Sky) — k — k^ and A^; < 0. The new 
effective Hamiltonian (5) has eigenvalues given by 



E+{5k) = ±hv3\ 



PD 



Pd_\ 
4 / 



4{Sky) 



(6) 

So, when -po/^ < < Q or Q < dy < pc/(4/3), there 
are still two Dirac cones at fcco + ^o and kni+Ai. When 
Xx reaches to a critical value of —pd/4:, the Hamiltonian 
is given by fiV3Tx{5kxY /p£, + 2hv3Ty6ky so that the dis- 
persion along kx direction is massive while one along ky 

direction is still massless. Two other cones at A:£i2 + A2 
and fc£)3 + A3 move away from each other in ±ky di- 
rection maintaining their anisotropic Dirac cone shapes. 
Therefore, it can be seen that the cones at kno and kui 
merge together when A = (— p_d/4, 0) [Fig. 4(c)]. The 
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corresponding sliding distance of dy is about 0.3% of Cc 
{dy = pd/{4:0) ~ 0.004A). The low energy bands already 
change their topology from the original structure under 
extremely small sliding. For further sliding, A^, < — P£)/4 
{dy > po/i'^li)), the merged cone eventually disappears 
and the spectrum of Eq. (6) develops an energy gap at 
k+ as shown in Fig. 4(d). The gap is linear with slid- 
ing distance as given by 2hv3\\x +pd/M — 2hv3l3dy ^ 
0.9 X [dy/ac] eV. The opening of energy gap signals a 
pair annihilation of two massless Dirac electrons with 
the opposite 'topological charges' of ±1 since the two 
Dirac cones at koo and koi have the Berry's phase of 
±7r respectively. This also confirms the Berry's phase 
conservation since the remaining two anisotropic Dirac 
cones at fc£)2 and fc/js give the total 4>b of 27r. 

Next, for the sliding along ~Si direction, the effec- 
tive Hamiltonian can be obtained by expanding Eq. (2) 
around fc_ = (pd/2,0): 



2Skx + 



{5k,f 



{6kyf 



Pd 



+ 



ipD 

4 



A. 



+2hvsT: 



{5kx){5ky) 



Pd 



(7) 



where 5k = (5kx,5ky) = k — k^ and A^; > 0. When 
< \x < 3pd/4, there are three Dirac cones at 
koi + Ai {i = 0,2,3). However, when A^; = 3pd/4, 
the three Dirac cones merge at fc_ and this Hamilto- 
nian gives an anomalous dispersion relation written by 



E^{5k) = ±hv:i^ p{5k) + 9'^{5k) where f{6k) = 26kx + 

{5kxf/pD - {5kyYlpD and g{5k) = 2{5kx){5ky) /po- 
This gives sickle-shaped constant energy contours which 
are consistent with our ah initio calculation results shown 
in Figs. 3(e) and (g). When A^: > 3pd/4, the effective 
Hamiltonian (7) does not develop any energy gap at all. 
Instead, when A^, ^ ipn /4, the effective Hamiltonian 
has a Dirac point at k'_ = {—pd/2 + \/\xPDt 0) and is 
given by 



V!_ ~ 2/IW3W —Tl5kx ± 2hv3 




<5ky, (8) 



where 5k = {5k,j.,5ky) ^ k — k'_. So, as sliding dis- 
tance increases along —5i, the anomalous Hamiltonian 
(7) gradually transforms to the anisotropic Dirac Hamil- 
tonian (8). This is quite contrary to the gapped spectrum 
(6) generated by sliding motion along 5i direction {T-L+ 
in Eq. (5)). The other cone at km + moves in —k^ 
direction maintaining its anisotropic dispersion relation. 
As Xx ^ 3p£)/4, the Dirac cone at koi + has an 
asymtotic shape as following, 



U'!_ ~ -2hvi 




-2hV3 




+ 1 r;5ky, (9) 



where 5k = (5kx,5ky) = k — kui —Ai. Therefore, in con- 
trast to the first case, no energy gap develops even when 
the sliding distance is increased further. Instead, the 
anomalous dispersion transforms to an anisotropic Dirac 
cone [Fig. 4(g)]. This phenomenon can be interpreted as 
a merging of two massless fermions of topological charge 
-f 1 with one of topological charge —1. As a result, a new 
fermion of topological charge +1 is generated. We note 
that the total (j)B of 27r is conserved since the new particle 
has (pB of TT. Hence, the topological charges of fermionic 
particles in BLG are strictly governed by Berry's phase 
conservation rule. 



VI. EFFECT OF SMALLEST INTERLAYER 
HOPPINGS 

When sliding occurs, the smallest interlayer interaction 
74 shown in Fig. 1(e) becomes anisotropic and depends 
on the pair distances between relevant carbon atoms in 
the top and bottom layer. This effect adds an additional 
Hamiltonian to Eq. (2), 



eflf 



71 
71 



k-{k + A)to 



{k + A)]7 



(k + X)], (10) 



where f4 = ^-^^ (See Appendix for derivation of Eq. 
(10)). Since the second term in Eq. (10) is twenty 
times smaller than the first term, we will neglect the sec- 
ond term hereafter. When sliding along ±5-1, the differ- 
ence between Dirac energies at kni and fc^o is given by 

27iii?ii4 



1(1 — 



1 



Without shding (A^^ = 0), the 



difference becomes Aeoh 



T^i : 



which indicates the 



hole and electron doped Dirac cone dX kui and fcuo re- 
spectively. This explains the hole and electron doped 
cones shown in our ah initio calculation resiilts [Fig. 
3(b)]. With sliding along 5\, the difference disappears 
when Xx < —po/^ as shown in Figs. 3(c) and 3(d) so 
that all Dirac cones are charge-neutral. Contrary to this, 
when sliding along ~5i, the difference changes its sign 
when Xx > Pd and decreases significantly as sliding dis- 
tance increases. So, the hole-doped Dirac cone at km 
changes to be electron-doped and the new Dirac cone to 
be hole-doped with increasing amount of doping as in- 
creasing sliding distance as shown in our first-principles 
calculations [Figs. 3 (c)-(h)]. We note that variations in 
74 do not affect any topological changes discussed so far. 



VII. DISCUSSION 

The direct signatures of the ETT can be readily mea- 
sured using various experiment methods. We showed 
that the tiny sliding lifts the degeneracy at Ep as the 
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1 1 

-(a) 


1 1 


^ \ 

1 1 


I 

1 1 




FIG. 5: First-principles calculations of total density of states 
of sliding BLG when ds = (a) 0, (b) 0.02(5i, (c) -0.02(5i, and 
(d) 0.012(53 — Si), respectively. The grey solid and dotted 
lines in (c) and (d) are the partial density of states projected 
onto Mix > i^s and < Ks region of the BZ in Fig. 3(e) and 
3(g), respectively, where component of the saddle 

point of each energy contour. 



number of Dirac cones is always reduced from four to 
two with significant increase of the saddle point energies 
and deformations of remaining Dirac cones. Therefore, 
first, high-resolution scanning tunneling microscopj*^^ 
can directly measure the changes in the density of states 
{D{E)). Between saddle point energies \E\ < Eg where 
= 7i^^3/(4t'o) without sliding, the total density of 
states per unit area is given by D{E) = \E\/{ttv'^). 
When BLG slides along ±Si, saddle point energies in- 
crease linearly as sliding distance increases, £s{dy) = 

fivsl^x—po/i] ~ hv3l3dy. When sliding along Si, the two 
Dirac cones disappear quickly and the other two remain- 
ing ones are anisotropic massless Dirac cones. So, the 
density of states in between saddle point energies ±es{dy) 
linearly depends on energy, D{E) ~ \E\. Our first- 
principles calculations indeed show drastic variations in 



the position of van Hove singularity (vHS) below and 
above the charge neutral point upon sliding (Fig. 5). In 
between the two vHSs, a linear [Figs. 5(a) and 5(b)] or a 
mostly square-root dependence of D{E) [Figs. 5(c) and 
(d)] appears as the sliding direction is changed, which is 
a unique feature of two-dimensional materials^^ii^. 

Second, the Landau level (LL) spectrum for a small 
perpendicular magnetic field (B) also exhibits distinctive 
dependence on the sliding direction. By using semiclas- 
sical quantization rule under the perpendicular magnetic 
fields (B), S{E) = ^^^(n + j) and dsSiE) = Att^D{E) 
where S{E) is an area of closed orbit of electron and 
7 — 1/2 — /2tv^^'^^^, we can immediately confirm 
that the Landau level (LL) spectrum under a small mag- 
netic field is given by En ^ ±(i?n)^/^ in the case of 

shding BLG along = 0.02(5i. When BLG slides by ei- 
ther — 0.02(5i or 0.012(53 — (5i), the anomalous Dirac cone 
shape results in D{E) ^ \fE being similar to the previous 
study on the density of states of semi-Dirac cone (mas- 
sive in one direction and massless in the other )^-'"^. By 
using the rule above and D{E) ^ \/E, the LL spectrum 
is given by En ±{Bn)'^/^. It is noticeable that, irre- 
spective of sliding direction, the zeroth LL exists at zero 
energy since the topological charge conservation (Berry's 
phase conservation) enforce the existence of at least one 
massless modes in the system. 

Third, an ideal minimal conductivity of 24e^/(7r^?i) at 
the charge neutral pointAi^i- will decrease quickly when 
sliding occurs (e is the electron charge). The conductiv- 
ity of anisotropic massless Dirac fermions with a disper- 
sion, VxTxkx + VyTyky IS glvcu hj (7 = whcrc 
5 is a degeneracy factor (5 = 4 if considering valley and 
spin degeneracies) and e is an electron charge. So, by 
using Eqs (8) and (9), the conductivity in wide BLG 
sliding along —Si {dy ^ pn) can be calculated easily, 

a = {I-pd/XxV^ ^ ^ (I+Pd/Xx). The conduc- 
tivity approaches 8e^/(7r^?i) as sliding distance increases. 

In summary, we show that the topology of energy 
bands of BLG changes significantly if sliding of extremely 
small distance occurs in any direction. The effective non- 
Abelian background gauge potential can be generated by 
sliding motion and is shown to play an important role 
in dictating the characteristics of sliding induced ETT. 
Hence, the ETT driven by non-Abelian gauge fields that 
are thought to be possible in cold atomic gas^i^i^ or 
similar effects in high energy physics^ can be realizable 
in sliding BLG. 

Note added in proof After submission, we became 
aware of related works on similar systems from other 
groups^"— 
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Appendix: Derivation of Eqs. (2) and (10) 

A single particle Hamiltonian of BLG in Bernal stack- 
ing, H = Hq + Hi + H2, can be written as 



Ho = 
Hi = 

H2 = 



-7i^a|(f)62(r) 



-73 E 4(^^+ 5j)h{f} + (h.c), (A.2) 



-74 



-74^6l(r-fl,>2(r) + (h.c), (A.3) 

where is a Hamiltonian for intra-layer interactions 
in each SLG and Hi and H2 are for inter-layer inter- 
actions between two SLG. Here a\{aa) and b\{ha) are 
the creation (annihilaton) operator of 7r-electron located 
at A- and i?-sublattice of layer a (= 1,2) respectively, 
r = raai + 110,2 (™ and n are integers), ai(= ^3 — 5i) 
and (22 (= ^2 — 5i) are unit vectors of hexagonal lattice 
of SLG^iiSrJ^ii^. 7g is the intra-layer nn hopping param- 
eter, and 71 and 73 are for the nn and next nn (nun) 
inter-layer hoppings parameters, respectively^iiSriSiii^. 

Using Fourier transformations of aa(r) — 
-l^J2^^e-'p-ra^p and 6„(r) = ^ .6"'^^-%^^ (N 

is a total number of unitcells), the total Hamilto- 
nian (H) in Eqs. (A.1)-(A.3) can be written in a 
matrix form such a.s H — X^p^^^p^p ^'^^ ^ 6.e\d of 
% = (bif, a2p, oijy, 62p)^ where 



Hp — 



Qip) Up) Up)' 

Up) aip) QiP) 

Coip) Up) -71 

va(p) Up) -71 . 



(A.4) 



and Up) = -7aEje'P'^ (" = 0,3,4, j = 1,2,3). 
When we expand Eq. (A.4) around if-point by using 
p = k + K (|fc| « \K\) and K = (^ ^^,0), 

Uk + K) ~la e'"^-'^ -ilo^Y^k- 5,e'^-'^ 
j j 
= + iky), (A. 5) 

where Va = ^laO-c (<^ = 0, 3, 4), S is the Planck constant 
and hk = h{kx,ky) is the crystal momentum from K- 



point. The total Hamiltonian near ivT-point for a field 
^j: — {bift, fljfci ^iki ^2fc)^ '■^^^^ written as 



/ W3fc_ Vok+ U4fc+^ 

v^k^ Vik^ Vfjk^ 

vok^ v^k+ —71 

\w4fc_ wofc+ —71 



(A.6) 



where k± — k^ i iky. The effective Hamiltonian (iJeff 
-ff^ff) on the low energy electronic structures for a field 
K 

k 



^I'r = (&i^,a2fe)'^' ^"^^ described by 



Hcff 



ft^^3 



k^ 
k+ 



2h'^VQV4 

71 

+ ^ 
71 







71 


k+k- 





fc2 



k\ 









(A.7) 



(A.8) 



Here we decompose the effective Hamiltonian into "Heff 
and "Hpg where the latter is quite small compared to the 
former. By using Pauli spin matrices, tq = (J^*)' — 
and Ty = ( ■ 7), the above Eqs. (A.7) and (A.8) 
can be written in compact forms. 



'H-cB 



nv3f-k+!^{T*-k)T^{T*-k), (A.9) 

71 

2h'^VoV4 I ^ 2 



M'ro 
71 

+ ^(^-fc)r.(r*.A;), 

71 



(A.IO) 



where r = (t^, Ty) and is its complex conjugate. Here- 
after, we will neglect the smallest hopping parameter 74 
(Eqs. (A.8) and (A.IO)) and discuss its role later. 

The effective Hamiltonian in Eq. (A.9) gives energy 
eigenvalues 



E{k) = ±hv:ik 



1 + — e^"t"' 

PD 



(A.ll) 



'^{ky/kx) and 



PD = 



where k — \k\, (j)k — tan 

7i^'3/(ft^'o) = 2/(3ac)7i73/7o- 

When the layer 2 (bottom layer) slides against the layer 

1 (top) along ds = {dx,dy) {\ds\ < Oc) [Figs. 1(a) and 
(b)], the constant nnn inter-layer interaction (73) now 
depends on the carbon pair distance in layer 1 and 2 [Fig. 
1(d)] and nn inter-layer interaction (71) decreases to 7^. 
Since we focus on an extremely small sliding distance, we 
assume variation of the nnn inter-layer interaction such 



as 73 lz{5i) ^ 736 



73(1 - I55i ■ ds) where i 



1, 2, 3 and /3 is a positive real constant. When we expand 
Eq. (A.4) around /T-point including sliding effect, intra- 
layer interactions remain the same as before. However, 
unlike the expansion procedure without sliding shown in 
Eq. (A. 5), the interlayer interaction with sliding can be 
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expanded up to a leading order of k and d as 



- ^3S^ ■ ds)il + ik ■ Sj)e' 



K-Si 



~ hv3{kx + iky) + hv3l3{id.j: - dy). (A. 12) 

Hence, the sliding vector ds plays a role of shifting the 
nnn intra-layer Hamiltonian in momentum space like a 
constant vector potential. With 74 neglected, the total 
Hamiltonian near ii'-point in Eq. (A. 6) now transforms 
to 



H-, = h 



I Vi{k- - A_) VQk^ 

i;3(fc+ — A+) wofc- 

wofc_ -7^ 

V wofc+ -7^ 

(A.13) 

where A± = T idx). Here we neglect an overall 

phase shift of e^"^''^ because of \ds\ ^ Oc. The effective 
low energy Hamiltonian of Eq. (A. 9) also changes to 

Heff ^ hv^T- (fc - A) + :^(f* . fc)r,(r' • fc) (A.14) 
7i 

where A = {\x-,\) = P{dy,~dx) = Pids x k^) {kz = 

Now, let us consider the smallest interlayer interac- 
tion 74 with sliding. When sliding occurs, 74 shown in 
Fig. 1(e) becomes anisotropic and depends on the pair 



distances between relevant carbon atoms in the top and 
bottom layer. We find that 74 — > ji{6i) ~ 74e''^' '^» ~ 
74(1 + /35i ■ ds) where i = 1,2,3. Here we assume the 
same coefficient /3 of 73 ((5i) for calculation convenience 
and note that it does not change the main conclusions of 
the paper. Like Eq. (A. 12) for 73, we can expand ^4 as 



Uk + K) ^ - Y^j,iS,)e^(^+^y'^ 
j 

- -5I'T4(l + /3<5. •rf;)(l + *fc-<5j)e" 

j 

~ hvi{kx + iky) — hv4f3{idx — dy) 



KSi 



hv/^k+ + hvi\+. 



(A.15) 



With 74 included, now the total Hamiltonian near K- 
point in Eq. (A.13) has an additional term. 



A:+ + A^ 

fc_ + A_ 

k+ + \+ 
vfc_ + A_ 



Then, Eq. (A. 10) changes to 



(A.16) 



^off 



71 



k-{k + X)to 



+ . (k + \)]T,[r ■ {k + X)]. (A.17) 

71 
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